package scu.maqiang.fes;

import java.sql.Array;
import java.util.*;
import java.util.stream.IntStream;

import scu.maqiang.mesh.*;
import scu.maqiang.numeric.Complex;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DoubleIntegerPair;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.Pair;
import scu.maqiang.numeric.ParamCheck;
import scu.maqiang.numeric.PointInElement;
import scu.maqiang.numeric.SCMatrix;
import scu.maqiang.numeric.SREntry;
import scu.maqiang.numeric.SRMatrix;
import scu.maqiang.numeric.Triple;

/**
 * 抽像有限元空间，为所有有限元空间的基类
 * 
 * @author 马强
 * @version 1.0
 * @since 1.0
 */
abstract class FES {

	/**
	 * 抽像有限元构造函数
	 * 
	 * @param mesh 网格类
	 */
	public FES(Mesh mesh) {
		this.mesh = mesh;
		nE = mesh.getNt();
		nB = mesh.getNb();
		nDim = mesh.getnDim();
		nPerEle = mesh.getnPerElement();
		nPerBoundary = mesh.getnPerBoundary();
		nBand = mesh.elementDiff() + 1;
		coordNode = new double[nDim];
		coordRegion = new double[nPerEle][nDim];
		coordBoundary = new double[nPerBoundary][nDim];
	}

	/**
	 * 构造函数
	 * 
	 * @param mesh       网格类
	 * @param dofPerNode 整数，每个结点自由度个数
	 */
	public FES(Mesh mesh, int dofPerNode) {
		this(mesh);
		this.dofPerNode = dofPerNode;
		ndof = dofPerNode * mesh.getNv();
		dofPerElement = dofPerNode * nPerEle;
		dofPerBoundary = dofPerNode * nPerBoundary;
		nBand = (mesh.elementDiff() + 1) * dofPerNode;
		initEleMatrixVector();
		int idx = 0;
		int nv = mesh.getNv();
		dofs = new int[mesh.getNv()][];
		for(int i = 0; i < nv; i++) {
			int[] dof = new int[dofPerNode];
			for(int j = 0; j < dofPerNode; j++) {
				dof[j] = idx++;
			}
			dofs[i] = dof;
		}
	}

	/**
	 * 得到有限元空间自由度个数
	 * 
	 * @return 自由度个数
	 */
	public int GetNdof() {
		return ndof;
	}

	/**
	 * 得到有限元空间网格
	 * 
	 * @return 网格对象
	 */
	public Mesh getMesh() {
		return mesh;
	}

	/**
	 * 得到有限元空间产生的稀疏矩阵半带宽
	 * 
	 * @return
	 */
	public int getnBand() {
		return nBand;
	}

	public int[] getDofFromDofs(int[] ele) {
		int totalNum = 0;
		int localLen = ele.length;
		for(int i = 0; i < localLen; i++) {
			totalNum += dofs[ele[i]].length;
		}
		int[] localDof = new int[totalNum];
		int idx = 0;
		for(int i = 0; i < localLen; i++) {
			int[] dof = dofs[ele[i]];
			for(int j = 0, nj = dof.length; j < nj; ++j) {
				localDof[idx++] = dof[j];
			}
		}
		return localDof;
	}

	public int[] getDofFromDofs(int[] ele, int[] range) {
		int localLen = ele.length;
		int rangeLen = range.length;
		int totalNum = ele.length * range.length;
		int[] localDof = new int[totalNum];
		int idx = 0;
		for(int i = 0; i < localLen; i++) {
			for(int j = 0; j < rangeLen; j++) {
				localDof[idx++] = dofs[ele[i]][range[j]];
			}
		}
		return localDof;
	}

	public int[] getDofByElement(int[] ele) {;
		ArrayList<Integer> dofInEle = new ArrayList<>();
		for(int idx: ele) {
			for(int localDof : dofs[idx]) {
				dofInEle.add(localDof);
			}
		}
		return dofInEle.stream().mapToInt(Integer::valueOf).toArray();
	}


	public int[] getDofFromNodes(int[] nodes) {
		ArrayList<Integer> result = new ArrayList<>();
		for(int i = 0, n = nodes.length; i < n; i++) {
			int[] dof = dofs[nodes[i]];
			for(int j = 0, nj = dof.length; j < nj; j++) {
				result.add(dof[j]);
			}
		}
		return result.stream().mapToInt(Integer::valueOf).toArray();
	}

	public ArrayList<Integer> getDofFromNodes(ArrayList<Integer> nodes) {
		ArrayList<Integer> result = new ArrayList<>();
		for(int i = 0, n = nodes.size(); i < n; i++) {
			int[] dof = dofs[nodes.get(i)];
			for(int j = 0, nj = dof.length; j < nj; j++) {
				result.add(dof[j]);
			}
		}
		return result;
	}

	public ArrayList<Integer> getDofFromNodes(HashSet<Integer> nodes) {
		ArrayList<Integer> result = new ArrayList<>();
		for(Integer node: nodes) {
			int[] dof = dofs[node];
			for(int j = 0, nj = dof.length; j < nj; j++) {
				result.add(dof[j]);
			}
		}
		return result;
	}

	public ArrayList<Integer> getDofFromNodes(HashSet<Integer> nodes, int[] range) {
		ArrayList<Integer> result = new ArrayList<>();
		for(Integer node: nodes) {
			int[] dof = dofs[node];
			for(int j = 0, nj = range.length; j < nj; j++) {
				result.add(dof[range[j]]);
			}
		}
		return result;
	}

	public ArrayList<Integer> getDofFromNodes(ArrayList<Integer> nodes, int[] range) {
		ArrayList<Integer> result = new ArrayList<>();
		for(int i = 0, n = nodes.size(); i < n; i++) {
			int[] dof = dofs[nodes.get(i)];
			for(int j = 0, nj = range.length; j < nj; j++) {
				result.add(dof[range[j]]);
			}
		}
		return result;
	}


	public int[] getDofsByNode(int idx) {
		return dofs[idx];
	}

	public int[] getDofsByNodes(int[] idxes) {
		ArrayList<Integer> dofInNodes = new ArrayList<>();
		for(int node: idxes) {
			for(int localDof : dofs[node]) {
				dofInNodes.add(localDof);
			}
		}
		return dofInNodes.stream().mapToInt(Integer::valueOf).toArray();
	}

	public int[] getDofsByNodes(int[] idxes, int[] range) {
		ArrayList<Integer> dofInNodes = new ArrayList<>();
		for(int node: idxes) {
			int[] localDof = dofs[node];
			for(int i = 0, n = range.length; i < n; i++) {
				dofInNodes.add(localDof[range[i]]);
			}
		}
		return dofInNodes.stream().mapToInt(Integer::valueOf).toArray();
	}

	public int[] getDofsByNodes(HashSet<Integer> idxes) {
		ArrayList<Integer> dofInNodes = new ArrayList<>();
		for(int node: idxes) {
			for(int localDof : dofs[node]) {
				dofInNodes.add(localDof);
			}
		}
		return dofInNodes.stream().mapToInt(Integer::valueOf).toArray();
	}

	public int[] getDofsByNodes(HashSet<Integer> idxes, int[] range) {
		ArrayList<Integer> dofInNodes = new ArrayList<>();
		for(int node: idxes) {
			int[] localDof = dofs[node];
			for(int i = 0, n = range.length; i < n; i++) {
				dofInNodes.add(localDof[range[i]]);
			}
		}
		return dofInNodes.stream().mapToInt(Integer::valueOf).toArray();
	}

	public int[] getDofsByElement(int idx) {
		return getDofsByNodes(mesh.getElement(idx));
	}

	public int[] getDofsByElement(int idx, int[] range) {
		return getDofsByNodes(mesh.getElement(idx), range);
	}

	public int[] getDofsByBoundary(int idx) {
		return getDofsByNodes(mesh.getBoundary(idx));
	}

	public int[] getDofsByBoundary(int idx, int[] range) {
		return getDofsByNodes(mesh.getBoundary(idx), range);
	}

	public void printDofs() {
		int n = mesh.getNv();
		int idx = 0;
		for(int[] dof : dofs) {
			int len = dof.length;
			System.out.print(idx + ": ");
			for(int i = 0; i < len; i++) {
				System.out.print(dof[i] + "  ");
			}
			System.out.println("");
			idx++;
		}
	}

	public int[][] getDofs() {
		int[][] result = new int[dofs.length][];
		for(int i = 0, n = dofs.length; i < n; i++) {
			result[i] = dofs[i].clone();
		}
		return result;
	}

	public void setDofsFromMeshNodeInfo() {
		int index = 0;
		ndof = 0;
		for(int i = 0, n = mesh.getNv(); i < n; ++i) {
			int nDof_local = mesh.getNodeLabel(i);
			int[] dof = new int[nDof_local];
			for(int j = 0; j < nDof_local; j++) {
				dof[j] = index++;
			}
			dofs[i] = dof;
			ndof += nDof_local;
		}
	}

	/**
	 * 在得到有限元解基础之上，计算区域某些点处的解的值，如果在网格结点，则返回网格结点值，如果不在结点处，
	 * 则使用插值方法返回值，同时返回该坐标点位于的单元结点编号数组以及在各结点上的插值系数数组
	 * 
	 * @param u      有限元空间标量函数解
	 * @param xyz    有限元空间点坐标数组
	 * @param eleIdx 坐标点xyz所位于单元的结点编号数组
	 * @param eleVal 坐标点xyz所位于单元的结点上插值系数
	 * @param pie    函数接口, 判断坐标是否在单元内
	 * @return 有限元解在该点处的函数值
	 */
	protected double getValue0(double[] u, double[] xyz, int[] eleIdx, double[] eleVal, PointInElement pie) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double[] uu = new double[nPerEle];
		double value = 0;
		for (int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			if (pie.action(coordRegion, xyz, eleVal)) {
				for (int j = 0; j < nPerEle; j++) {
					uu[j] = u[ele[j]];
					eleIdx[j] = ele[j];
					value += uu[j] * eleVal[j];
				}
				return value;
			}
		}
		// throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + "
		// ) is not in the domain!");
		return 1.0e30;
	}
	
	/**
	 * 二分查找得到xyz位于哪一个单元之内
	 * @param u
	 * @param xyz
	 * @param eleIdx
	 * @param eleVal
	 * @param pie
	 * @return
	 */
	protected double getValue1(double[] u, double[] xyz, int[] eleIdx, double[] eleVal, PointInElement pie) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double[] uu = new double[nPerEle];
		double value = 0;
		for (int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			if (pie.action(coordRegion, xyz, eleVal)) {
				for (int j = 0; j < nPerEle; j++) {
					uu[j] = u[ele[j]];
					eleIdx[j] = ele[j];
					value += uu[j] * eleVal[j];
				}
				return value;
			}
		}
		// throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + "
		// ) is not in the domain!");
		return 1.0e30;
	}
	
	protected double getValue0(double[] u, double[] xyz, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double[] uu = new double[nPerEle];
		double value = 0;
		PointInElement pie = mesh.getPointJudgeFunc();
		for (int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			if (pie.action(coordRegion, xyz, eleVal)) {
				for (int j = 0; j < nPerEle; j++) {
					uu[j] = u[ele[j]];
					eleIdx[j] = ele[j];
					value += uu[j] * eleVal[j];
				}
				return value;
			}
		}
		// throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + "
		// ) is not in the domain!");
		return 1.0e30;
	}
	
	/**
	 * 在得到有限元解基础之上，计算区域某些点处的解的值，如果在网格结点，则返回网格结点值，如果不在结点处， 则使用插值方法返回值
	 * 
	 * @param u   有限元空间标量函数解
	 * @param xyz 有限元空间点坐标数组
	 * @param pie 函数接口, 判断坐标是否在单元内
	 * @return 有限元解在该点处的函数值
	 */
	protected double getValue0(double[] u, double[] xyz, PointInElement pie) {
		ParamCheck.checkEqual(u.length, this.ndof);
		int[] eleIdx = new int[nPerEle];
		double[] eleVal = new double[nPerEle];
		return getValue0(u, xyz, eleIdx, eleVal, pie);
	}

	/**
	 * 在得到有限元解基础之上，计算区域某些点处的解的值，如果在网格结点，则返回网格结点值，如果不在结点处，
	 * 则使用插值方法返回值，同时返回该点位于的单元编号，并且返回该单元的结点编号数组以及在各结点上的插值系数数组
	 * 
	 * @param u      有限元空间标量函数解
	 * @param xyz    有限元空间点坐标数组
	 * @param data   实数整数对对象，对象中第一个分量为函数值，第二个分量返回单元编号
	 * @param eleIdx 坐标点xyz所位于单元的结点编号数组
	 * @param eleVal 坐标点xyz所位于单元的结点上插值系数
	 * @param pie    函数接口, 判断坐标是否在单元内
	 */
	protected void getValue0(double[] u, double[] xyz, DoubleIntegerPair data, int[] eleIdx, double[] eleVal,
			PointInElement pie) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double value = 0;
		int[] ele = mesh.getElement(data.second);
		mesh.getCoordPerElement(ele, coordRegion);
		if (pie.action(coordRegion, xyz, eleVal)) {
			for (int j = 0; j < nPerEle; j++) {
				// eleIdx[j] = ele.idx[j];
				value += u[ele[j]] * eleVal[j];
			}
			data.first = value;
			return;
		} else {
			for (int i = 0, n = mesh.getNt(); i < n; i++) {
				ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				if (pie.action(coordRegion, xyz, eleVal)) {
					for (int j = 0; j < nPerEle; j++) {
						eleIdx[j] = ele[j];
						value += u[ele[j]] * eleVal[j];
					}
					data.first = value;
					data.second = i;
					return;
				}
			}

		}
		throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
	}
	
	protected void getValue0_V1(double[] u, double[] xyz, DoubleIntegerPair data, int[] eleIdx, double[] eleVal,
			PointInElement pie) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double value = 0;
		int[] ele = mesh.getElement(data.second);
		mesh.getCoordPerElement(ele, coordRegion);
		if (pie.action(coordRegion, xyz, eleVal)) {
			for (int j = 0; j < nPerEle; j++) {
				// eleIdx[j] = ele.idx[j];
				value += u[ele[j]] * eleVal[j];
			}
			data.first = value;
			return;
		} else {
			Triple<Integer, int[], double[]> info = mesh.pointInMesh(xyz, pie);
			ele = info.getSecond();
			eleVal = info.getThird();
			for (int j = 0; j < nPerEle; j++) {
				eleIdx[j] = ele[j];
				value += u[ele[j]] * eleVal[j];
			}
			data.first = value;
			data.second = info.getFirst();
		}
		throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
	}

	/**
	 * 提取向量值函数的各分量解
	 * 
	 * @param solu 问题解向量
	 * @param u    二维数组，行数为有限元网格上每个结点的自由度个数，列数为有限元网格结点数
	 */
	public void extractComponent(double[] solu, double[][] u) {
		// ParamCheck.checkEqual(dofPerNode, u.length);
		int m = u[0].length, idx;
		ParamCheck.checkEqual(solu.length, m * dofPerNode);
		for (int i = 0; i < m; i++) {
			idx = i * dofPerNode;
			for (int j = 0; j < dofPerNode; j++) {
				u[j][i] = solu[idx + j];
			}
		}
	}

	/**
	 * 初始化单元矩阵与向量，该方法在子类中调用，以便在整体刚度矩阵的组装过程中提供通用的单元矩阵与单元向量
	 */
	protected void initEleMatrixVector() {
		eleMatrixRegion = new double[dofPerElement][dofPerElement];
		eleMatrixBoundary = new double[dofPerBoundary][dofPerBoundary];
		eleVecRegion = new double[dofPerElement];
		eleVecBoundary = new double[dofPerBoundary];
	}

	/**
	 * 在整体向量中组装单元向量
	 * 
	 * @param ele    单元向量
	 * @param dof    自由度，需与单元向量长度相同
	 * @param global 整体向量
	 */
	void assembleVector(double[] ele, int[] dof, double[] global) {
		int ne = ele.length;
		ParamCheck.checkEqual(ne, dof.length);
		for (int i = 0; i < ne; i++) {
			global[dof[i]] += ele[i];
		}
	}

	/**
	 * 向整体复值型向量中组装单元复向量
	 * @param ele 单元向量
	 * @param dof 自由度, 需与单元向量长度相同
	 * @param global 整体向量
	 */
	void assembleCVector(Complex[] ele, int[] dof, Complex[] global) {
		int ne = ele.length;
		ParamCheck.checkEqual(ne, dof.length);
		for (int i = 0; i < ne; i++) {
			global[dof[i]].real += ele[i].real;
			global[dof[i]].image += ele[i].image;
		}
	}
	

	/**
	 * 得到指定边界所有结点相关的组装得到的稀疏矩阵对角线元素
	 * 
	 * @param A      稀疏矩阵
	 * @param direct 方向
	 * @param label  网格边界编号
	 * @return 得到实数组与整数组对, 实数组存储稀疏矩阵对角线元素, 整数组存储对角线坐标
	 */
	public Pair<double[], int[]> getDiagElementsFromBoundaries(SRMatrix A, Direct direct, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		int offset = direct.getValue();
		int idx = 0;
		double[] result;
		int[] dof;
		if (offset >= 0) {
			result = new double[bcSet.size()];
			dof = new int[bcSet.size()];
			for (int ele : bcSet) {
				dof[idx] = ele * dofPerNode + offset;
				result[idx] = A.getElement(dof[idx], dof[idx]);
				idx++;
			}
		} else {
			result = new double[bcSet.size() * dofPerNode];
			dof = new int[bcSet.size() * dofPerNode];
			for (int ele : bcSet) {
				for (int j = 0; j < dofPerNode; j++) {
					dof[idx] = ele * dofPerNode + j;
					result[idx] = A.getElement(dof[idx], dof[idx]);
					idx++;
				}
			}
		}
		return new Pair<>(result, dof);
	}

	/**
	 * 施加周期边界条件
	 * @param A 整体矩阵
	 * @param RHS 右端项
	 * @param idxPair 指标对, 行数为边界自由度数, 列数为2
	 */
	public void applyBC_Periodic(SRMatrix A, double[] RHS, int[][] idxPair) {
		for(int[] pair:idxPair) {
			int[] dofPair = this.getDof(pair);
			for(int k = 0; k < dofPerNode; k++) {
				int dof1 = dofPair[k];
				int dof2 = dofPair[k + dofPerNode];
				RHS[dof2] += RHS[dof1];
				RHS[dof1] = 0.0;
				for(SREntry ele: A.mat[dof1]) {
					A.addElement(dof2, ele.col, ele.val);
				}
				for(int i = 0, n = A.getM(); i < n; i++) {
					double temp = A.getElement(i, dof1);
					A.removeElement(i, dof1);
					A.addElement(i, dof2, temp);
				}
				//System.out.println(A.isSymmetric());
				//A.removeSymmetricIJ(dof1);
				A.removeRow(dof1);
				//System.out.println(A.isSymmetric());
				A.setElement(dof1, dof1, 1.0);
				A.setElement(dof1, dof2, -1.0);
				A.setElement(dof2, dof1, -1.0);
				A.addElement(dof2, dof2, 1.0);
				//System.out.println(A.display());
			}
		}
		//System.exit(1);
	}
	
	public void applyBC_Periodic(SRMatrix A, int[][] idxPair) {
		for(int[] pair:idxPair) {
			int[] dofPair = this.getDof(pair);
			for(int k = 0; k < dofPerNode; k++) {
				int dof1 = dofPair[k];
				int dof2 = dofPair[k + dofPerNode];
				for(SREntry ele: A.mat[dof1]) {
					A.addElement(dof2, ele.col, ele.val);
				}
				for(int i = 0, n = A.getM(); i < n; i++) {
					double temp = A.getElement(i, dof1);
					A.removeElement(i, dof1);
					A.addElement(i, dof2, temp);
				}
				//System.out.println(A.isSymmetric());
				//A.removeSymmetricIJ(dof1);
				A.removeRow(dof1);
				//System.out.println(A.isSymmetric());
				A.setElement(dof1, dof1, 1.0);
				A.setElement(dof1, dof2, -1.0);
				A.setElement(dof2, dof1, -1.0);
				A.addElement(dof2, dof2, 1.0);
				//System.out.println(A.display());
			}
		}
	}
	
	public void applyBC_Periodic(double[] RHS, int[][] idxPair) {
		for(int[] pair:idxPair) {
			int[] dofPair = this.getDof(pair);
			for(int k = 0; k < dofPerNode; k++) {
				int dof1 = dofPair[k];
				int dof2 = dofPair[k + dofPerNode];
				RHS[dof2] += RHS[dof1];
				RHS[dof1] = 0.0;
			}
		}
	}
	
	/**
	 * 施加周期性边界条件
	 * @param A
	 * @param RHS
	 * @param direct
	 * @param label1
	 * @param label2
	 */
	public void applyBC_Periodic(SRMatrix A, double[] RHS, Direct direct, int label1, int label2) {
		int[][] idxPair = mesh.extractPeriodicBoundary(label1, label2, direct);
//		System.out.println(MVO.toString(idxPair));
		applyBC_Periodic(A, RHS, idxPair);
	}

	public void applyBC_Periodic(SRMatrix A, Direct direct, int label1, int label2) {
		int[][] idxPair = mesh.extractPeriodicBoundary(label1, label2, direct);
		applyBC_Periodic(A, idxPair);
	}
	
	public void applyBC_Periodic(double[] RHS, Direct direct, int label1, int label2) {
		int[][] idxPair = mesh.extractPeriodicBoundary(label1, label2, direct);
		applyBC_Periodic(RHS, idxPair);
	}
	
	public void applyBC_Periodic(SRMatrix A, double[] RHS, double[] tangent, int label1, int label2) {
		int[][] idxPair = mesh.extractPeriodicBoundary(label1, label2, tangent);
		applyBC_Periodic(A, RHS, idxPair);
	}

	public void applyBC_Periodic(SRMatrix A, Direct direct, double[] tangent, int label1, int label2) {
		int[][] idxPair = mesh.extractPeriodicBoundary(label1, label2, tangent);
		applyBC_Periodic(A, idxPair);
	}
	
	public void applyBC_Periodic(double[] RHS, Direct direct, double[] tangent, int label1, int label2) {
		int[][] idxPair = mesh.extractPeriodicBoundary(label1, label2, tangent);
		applyBC_Periodic(RHS, idxPair);
	}
	
	
	/**
	 * 使用乘大数法对线性方程组Ax = b 施加边界条件
	 * 
	 * @param A
	 * @param RHS
	 * @param direct
	 * @param val
	 * @param label
	 */
	public void applyBC_MBN(SRMatrix A, double[] RHS, Direct direct, double val, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		applyBC_MBN(A, RHS, direct, val, bcSet);
	}

	public void applyBC_MBN(SRMatrix A, double[] RHS, Direct direct, double val, HashSet<Integer> bcSet) {
		int offset = direct.getValue();
		if (offset >= 0) {
			bcSet.forEach(dof -> {
				int localIdx = dof * dofPerNode + offset;
				A.setElement(localIdx, localIdx, 1.0e30);
				RHS[dof * dofPerNode + offset] = val * 1.0e30;
			});
		} else {
			bcSet.forEach(dof -> {
				for (int j = 0; j < dofPerNode; j++) {
					int localIdx = dof * dofPerNode + j;
					A.setElement(localIdx, localIdx, 1.0e30);
					RHS[dof * dofPerNode + j] = val * 1.0e30;
				}
			});
		}
	}

	public void applyBC_MBN(SRMatrix A, double[] RHS, Direct direct, double val, HashSet<Integer> bcSet,
			double[] diagValue) {
		int offset = direct.getValue();
		if (offset >= 0) {
			ParamCheck.checkEqual(diagValue.length, bcSet.size());
			bcSet.forEach(dof -> {
				int localIdx = dof * dofPerNode + offset;
				A.setElement(localIdx, localIdx, 1.0e30);
				RHS[dof * dofPerNode + offset] = val * 1.0e30;
			});
		} else {
			ParamCheck.checkEqual(diagValue.length, bcSet.size() * dofPerNode);
			bcSet.forEach(dof -> {
				for (int j = 0; j < dofPerNode; j++) {
					int localIdx = dof * dofPerNode + j;
					A.setElement(localIdx, localIdx, 1.0e30);
					RHS[dof * dofPerNode + j] = val * 1.0e30;
				}
			});
		}
	}

	/**
	 * 通过乘大数法在特定网格边界上设置稀疏矩阵边界条件
	 * 
	 * @param A      整体稀疏矩阵
	 * @param direct 设置边界条件的自由度分量，如果问题为标量问题，可取direct值为X或ALL, 若问题为向量问题，则按实际情况取值
	 * @param label  网格边界标识号
	 */
	public void applyBC_MBN(SRMatrix A, Direct direct, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		int offset = direct.getValue();
		if (offset >= 0) {
			bcSet.forEach(dof -> {
				int localIdx = dof * dofPerNode + offset;
				A.setElement(localIdx, localIdx, 1.0e30);
			});
		} else {
			bcSet.forEach(dof -> {
				for (int j = 0; j < dofPerNode; j++) {
					int localIdx = dof * dofPerNode + j;
					A.setElement(localIdx, localIdx, 1.0e30);
				}
			});
		}
	}
	
	/**
	 * 通过乘大数法在特定网格边界上设置稀疏矩阵边界条件
	 * 
	 * @param A      整体稀疏矩阵
	 * @param direct 设置边界条件的自由度分量，如果问题为标量问题，可取direct值为X或ALL, 若问题为向量问题，则按实际情况取值
	 * @param label  网格边界标识号
	 */
	public void applyBC_MBN(SCMatrix A, Direct direct, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		int offset = direct.getValue();
		if (offset >= 0) {
			bcSet.forEach(dof -> {
				int localIdx = dof * dofPerNode + offset;
				A.setElement(localIdx, localIdx, new Complex(1.0e30, 0));
			});
		} else {
			bcSet.forEach(dof -> {
				for (int j = 0; j < dofPerNode; j++) {
					int localIdx = dof * dofPerNode + j;
					A.setElement(localIdx, localIdx, new Complex(1.0e30, 0));
				}
			});
		}
	}

	public void applyBC_MBN(SRMatrix A, double[] RHS, int[] dofs, double[] vals) {
		for(int i = 0, n = vals.length; i < n; i++) {
			A.setElement(dofs[i], dofs[i], 1.0e30);
			RHS[dofs[i]] = vals[i] * 1.0e30;
		}
	}

	public void applyBC_MBN(SRMatrix A, int[] dofs) {
		IntStream.of(dofs).forEach(dof -> A.setElement(dof, dof, 1.0e30));
	}

	public void applyBC_MBN(double[] RHS, int[] dofs, double[] vals) {
		for(int i = 0, n = vals.length; i < n; i++) {
			RHS[dofs[i]] = vals[i] * 1.0e30;
		}
	}

	public void applyBC_MBN(SCMatrix A, Complex[] RHS, int[] dofs, Complex[] vals) {
		for(int i = 0, n = vals.length; i < n; i++) {
			A.setElement(dofs[i], dofs[i], new Complex(1.0e30, 0.0));
			RHS[dofs[i]].real = vals[i].real * 1.0e30;
			RHS[dofs[i]].image = vals[i].image * 1.0e30;
		}
	}

	public void applyBC_MBN(SCMatrix A, int[] dofs) {
		IntStream.of(dofs).forEach(dof -> A.setElement(dof, dof, new Complex(1.0e30, 0)));
	}

	public void applyBC_MBN(Complex[] RHS, int[] dofs, Complex[] vals) {
		for(int i = 0, n = vals.length; i < n; i++) {
			RHS[dofs[i]].real = vals[i].real * 1.0e30;
			RHS[dofs[i]].image = vals[i].image * 1.0e30;
		}
	}


	/**
	 * 通过乘大数法设置在特定网格边界上右端项边界条件
	 * 
	 * @param RHS    代数方程组右端项
	 * @param direct 设置边界条件的自由度分量，如果问题为标量问题，可取direct值为X或ALL, 若问题为向量问题，则按实际情况取值
	 * @param val    要设置的边界值
	 * @param label  网格边界标识号
	 */
	public void applyBC_MBN(double[] RHS, Direct direct, double val, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		int offset = direct.getValue();
		if (offset >= 0) {
			bcSet.forEach(dof -> {
				RHS[dof * dofPerNode + offset] = val * 1.0e30;
			});
		} else {
			bcSet.forEach(dof -> {
				for (int j = 0; j < dofPerNode; j++) {
					RHS[dof * dofPerNode + j] = val * 1.0e30;
				}
			});
		}
	}

	public void applyBC_MBN(Complex[] RHS, Direct direct, Complex val, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		int offset = direct.getValue();
		if (offset >= 0) {
			bcSet.forEach(dof -> {
				RHS[dof * dofPerNode + offset].real = val.real * 1.0e30;
				RHS[dof * dofPerNode + offset].image = val.image * 1.0e30;
			});
		} else {
			bcSet.forEach(dof -> {
				for (int j = 0; j < dofPerNode; j++) {
					RHS[dof * dofPerNode + j].real = val.real * 1.0e30;
					RHS[dof * dofPerNode + j].image = val.image * 1.0e30;
				}
			});
		}
	}

	public void applyBC_MBN(Complex[] RHS, Direct direct, ScalarCFunc[] func, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		double[] xy = new double[this.nDim];
		int offset = direct.getValue();
		if (offset >= 0) {
			bcSet.forEach(dof -> {
				mesh.getCoordPerNode(dof, xy);
				Complex val = func[0].action(xy, 0, null);
				RHS[dof * dofPerNode + offset].real = val.real * 1.0e30;
				RHS[dof * dofPerNode + offset].image = val.image * 1.0e30;
			});
		} else {
			bcSet.forEach(dof -> {
				mesh.getCoordPerNode(dof, xy);
				for (int j = 0; j < dofPerNode; j++) {
					Complex val = func[j].action(xy, 0, null);
					RHS[dof * dofPerNode + j].real = val.real * 1.0e30;
					RHS[dof * dofPerNode + j].image = val.image * 1.0e30;
				}
			});
		}
	}

	public void applyBC_MBN_Disp(double[] RHS, Direct direct, double[] val, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		int offset = direct.getValue();
		if (offset > 0) {
			// ParamCheck.checkEqual(bcSet.size(), val.length);
			for (int dof : bcSet) {
				RHS[dof * dofPerNode + offset] = val[0] * 1.0e30;
			}
		} else {
			// ParamCheck.checkEqual(bcSet.size() * dofPerNode, val.length);
			for (int dof : bcSet) {
				for (int j = 0; j < dofPerNode; j++) {
					RHS[dof * dofPerNode + j] = val[j] * 1.0e30;
				}
			}
		}
	}

	public void applyBC_MBN(SRMatrix A, double[] RHS, int[] idxSet, double[] val, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		for (int i = 0, n = idxSet.length; i < n; i++) {
			int offset = idxSet[i];
			ParamCheck.checkBounds(offset, 0, dofPerNode);
			for (int dof : bcSet) {
				int localIdx = dof * dofPerNode + offset;
				A.setElement(localIdx, localIdx, 1.0e30);
				RHS[localIdx] = val[i] * 1.0e30;
			}

		}
	}

	public void applyBC_MBN(SRMatrix A, int[] idxSet, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		for (int i = 0, n = idxSet.length; i < n; i++) {
			int offset = idxSet[i];
			ParamCheck.checkBounds(offset, 0, dofPerNode);
			for (int dof : bcSet) {
				int localIdx = dof * dofPerNode + offset;
				A.setElement(localIdx, localIdx, 1.0e30);
			}
		}
	}

	public void applyBC_MBN(double[] RHS, int[] idxSet, double[] val, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		for (int i = 0, n = idxSet.length; i < n; i++) {
			int offset = idxSet[i];
			ParamCheck.checkBounds(offset, 0, dofPerNode);
			for (int dof : bcSet) {
				int localIdx = dof * dofPerNode + offset;
				RHS[localIdx] = val[i] * 1.0e30;
			}

		}
	}

	/**
	 * ͨ���˴������������ض�����߽����Ҷ���߽�����
	 * 
	 * @param RHS    ��װ�õ��Ĵ�����������Ҷ���
	 * @param direct ���ñ߽����������ɶȷ������������Ϊ�������⣬��ȡdirectֵΪX��ALL,
	 *               ������Ϊ�������⣬��ʵ�����ȡֵ
	 * @param val    Ҫ���õı߽�ֵ�������������ȱ�����Ҫ��������߽��ϵĽ�������ͬ
	 * @param label  ����߽��ʶ��
	 */
	public void applyBC_MBN(double[] RHS, Direct direct, double[] val, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		ParamCheck.checkEqual(bcSet.size(), val.length);
		int offset = direct.getValue();
		int i = 0;
		if (offset > 0) {
			ParamCheck.checkEqual(bcSet.size(), val.length);
			for (int dof : bcSet) {
				RHS[dof * dofPerNode + offset] = val[i++] * 1.0e30;
			}
		} else {
			ParamCheck.checkEqual(bcSet.size() * dofPerNode, val.length);
			for (int dof : bcSet) {
				for (int j = 0; j < dofPerNode; j++) {
					RHS[dof * dofPerNode + j] = val[i++] * 1.0e30;
				}
			}
		}
	}

	/**
	 * 通过乘大数法设置在特定网格边界上右端项边界条件，使用函数进行设置，
	 * 
	 * @param RHS    组装得到的代数方程组的右端项
	 * @param direct 设置边界条件的自由度分量，如果问题为标量问题，可取direct值为X或ALL, 若问题为向量问题，则按实际情况取值
	 * @param rf     区域函数，参数为边界结点的坐标值与参数param, 设置特定边界条件
	 * @param param  区域函数所带实型参数
	 * @param label  可变数组，网格边界标识号
	 */
	public void applyBC_MBN(double[] RHS, Direct direct, ScalarFunc[] rf, double[] param, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		int offset = direct.getValue();
		if (offset >= 0) {
			for (int dof : bcSet) {
				mesh.getCoordPerNode(dof, coordNode);
				RHS[dof * dofPerNode + offset] = rf[0].action(coordNode, 0, param) * 1.0e30;
			}
		} else {
			for (int dof : bcSet) {
				mesh.getCoordPerNode(dof, coordNode);
				for (int j = 0; j < dofPerNode; j++) {
					RHS[dof * dofPerNode + j] = rf[j].action(coordNode, 0, param) * 1.0e30;
				}
			}
			;
		}
	}
	
	public void applyBC_MBN(Complex[] RHS, Direct direct, ScalarCFunc[] rf, Complex[] param, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		int offset = direct.getValue();
		if (offset >= 0) {
			for (int dof : bcSet) {
				mesh.getCoordPerNode(dof, coordNode);
				RHS[dof * dofPerNode + offset].real = (rf[0].action(coordNode, 0, param)).real * 1.0e30;
				RHS[dof * dofPerNode + offset].image = (rf[0].action(coordNode, 0, param)).image * 1.0e30;
			}
		} else {
			for (int dof : bcSet) {
				mesh.getCoordPerNode(dof, coordNode);
				for (int j = 0; j < dofPerNode; j++) {
					RHS[dof * dofPerNode + j].real = (rf[j].action(coordNode, 0, param)).real * 1.0e30;
					RHS[dof * dofPerNode + j].image = (rf[j].action(coordNode, 0, param)).image * 1.0e30;
				}
			}
			;
		}
	}

	public void applyBC_MBN_Vector(double[] RHS, Direct direct, double[] val, int... label) {
		ParamCheck.checkEqual(val.length, mesh.getNv() * this.dofPerNode);
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		int offset = direct.getValue();
		if (offset > 0) {
			for (int dof : bcSet) {
				int offIdx = dof * dofPerNode + offset;
				RHS[offIdx] = val[offIdx] * 1.0e30;
			}
		} else {
			for (int dof : bcSet) {
				for (int j = 0; j < dofPerNode; j++) {
					int offIdx = dof * dofPerNode + j;
					RHS[offIdx] = val[offIdx] * 1.0e30;
				}
			}
		}
	}

	/**
	 * 对稀疏矩阵施加0边界条件
	 * 
	 * @param A      稀疏矩阵
	 * @param direct ���ñ߽����������ɶȷ������������Ϊ�������⣬��ȡdirectֵΪX��ALL,
	 *               ������Ϊ�������⣬��ʵ�����ȡֵ
	 * @param tp     �߽��������ͣ����ΪMBN����ʹ�ó˴����������ΪRCE����Ϊ������ȥ��
	 * @param label  网格边界编号
	 */
	public void applyBCZero(SRMatrix A, Direct direct, BCType tp, int... label) {
		if (tp == BCType.MBN) {
			applyBC_MBN(A, direct, label);
		} else if (tp == BCType.RCE) {
			HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
			int offset = direct.getValue();
			if (offset > 0) {
				for (int dof : bcSet) {
					int localIdx = dof * dofPerNode + offset;
					A.removeSymmetricIJ(localIdx);
					// A.removeRow(localIdx);
					A.setElementFast(localIdx, localIdx, 1.0);
				}
			} else {
				for (int dof : bcSet) {
					for (int j = 0; j < dofPerNode; j++) {
						int localIdx = dof * dofPerNode + j;
						A.removeSymmetricIJ(localIdx);
						// A.removeRow(localIdx);
						A.setElementFast(localIdx, localIdx, 1.0);
					}
				}
			}
		}

	}

	/**
	 * �Դ����������Ҷ�������0�߽�����
	 * 
	 * @param RHS    �����������Ҷ���
	 * @param direct ���ñ߽����������ɶȷ������������Ϊ�������⣬��ȡdirectֵΪX��ALL,
	 *               ������Ϊ�������⣬��ʵ�����ȡֵ
	 * @param label  ����߽��ʶ��
	 */
	public void applyBCZero(double[] RHS, Direct direct, int... label) {
		applyBC_MBN(RHS, direct, 0.0, label);
	}

	/**
	 * ��ϡ�����������������Ҷ���ͬʱ����0�߽�����
	 * 
	 * @param A      稀疏矩阵
	 * @param RHS    右端项
	 * @param direct ���ñ߽����������ɶȷ������������Ϊ�������⣬��ȡdirectֵΪX��ALL,
	 *               ������Ϊ�������⣬��ʵ�����ȡֵ
	 * @param tp     问题类型
	 * @param label  ����߽��ʶ��
	 */
	public void applyBCZero(SRMatrix A, double[] RHS, Direct direct, BCType tp, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		int offset = direct.getValue();
		if (tp == BCType.MBN) {
			if (offset >= 0) {
				bcSet.forEach(dof -> {
					int localIdx = dof * dofPerNode + offset;
					A.setElement(localIdx, localIdx, 1.0e30);
					RHS[localIdx] = 0.0;
				});
			} else {
				for (int dof : bcSet) {
					for (int j = 0; j < dofPerNode; j++) {
						int localIdx = dof * dofPerNode + j;
						A.setElement(localIdx, localIdx, 1.0e30);
						RHS[localIdx] = 0.0;
					}
				}
			}
		} else if (tp == BCType.RCE) {
			if (offset >= 0) {
				for (int dof : bcSet) {
					int localIdx = dof * dofPerNode + offset;
					//A.removeSymmetricIJ(localIdx);
					A.removeRowCol(localIdx);
					A.setElementFast(localIdx, localIdx, 1.0);
					RHS[localIdx] = 0.0;
				}
			} else {
				for (int dof : bcSet) {
					for (int j = 0; j < dofPerNode; j++) {
						int localIdx = dof * dofPerNode + j;
						//A.removeSymmetricIJ(localIdx);
						A.removeRowCol(localIdx);
						A.setElementFast(localIdx, localIdx, 1.0);
						RHS[localIdx] = 0.0;
					}
				}
			}
		}

	}

	/**
	 * 按边界自由度编号施加边界条件
	 * 
	 * @param A     稀疏矩阵
	 * @param RHS   右端项向量
	 * @param BCDOF 自由度数组
	 * @param val   施加的边界条件的值
	 * @param tp    边界条件施加类型
	 */
	public void applyBC_BCDOF(SRMatrix A, double[] RHS, int[] BCDOF, double[] val, BCType tp) {
		int n = BCDOF.length;
		int dof, colIdx;
		if (tp == BCType.MBN) {
			for (int i = 0; i < n; i++) {
				dof = BCDOF[i];
				RHS[dof] = val[i] * 1.0e30;
				A.setElement(dof, dof, 1.0e30);
			}
		} else if (tp == BCType.RCE) {
			for (int i = 0; i < n; i++) {
				dof = BCDOF[i];
				int len = A.mat[dof].size();
				for (int k = 0; k < len; k++) {
					colIdx = A.mat[dof].get(k).col;
					if (colIdx != i) {
						RHS[colIdx] -= val[i] * A.getElement(colIdx, dof);
						A.removeElement(colIdx, dof);
					}
				}
				A.removeRow(dof);
				A.setElementFast(dof, dof, 1.0);
				RHS[dof] = val[i];
			}
		}
	}
	
	
	public void applyBC_RCE(SRMatrix A, double[] RHS, Direct direct, double val, int... label) {
		HashSet<Integer> bcSet = mesh.extractBoundaryNodes(label);
		int offset = direct.getValue();
		SRMatrix AT = A.transpose();
		if(offset >= 0) {
			for(int ele : bcSet) {
				int dof = ele * dofPerNode + offset;
				double valII =  A.getElement(dof, dof);
				for(SREntry sele : AT.mat[dof]) {
					RHS[sele.col] -= valII * val;
				}
			}
		} else {
			
		}
	}
	
	/**
	 * 从函数中得到系数
	 * 
	 * @param coefFunc 标量函数数组
	 * @param coord    单元坐标数组ֵ
	 * @param label    单元材料标识
	 * @param param    标量函数参数
	 * @param coef     得到的系数数组
	 */
	public void getCoefFromFunc(ScalarFunc[] coefFunc, double[][] coord, int label, double[] param, double[][] coef) {
		for (int i = 0, ni = coefFunc.length; i < ni; i++) {
			for (int j = 0, nj = coef[i].length; j < nj; j++) {
				coef[i][j] = coefFunc[i].action(coord[j], label, param);
			}
		}
	}
	
	public void getCoefFromFunc(ScalarCFunc[] coefFunc, double[][] coord, int label, Complex[] param, Complex[][] coef) {
		for (int i = 0, ni = coefFunc.length; i < ni; i++) {
			for (int j = 0, nj = coef[i].length; j < nj; j++) {
				coef[i][j] = coefFunc[i].action(coord[j], label, param);
			}
		}
	}
	
	/**
	 * 得到某个结点的自由度
	 * 
	 * @param pt     网格结点编号
	 * @param direct 方向
	 * @return 得到结点在指定direct方向下的自由度
	 */
	public int[] getDof(int pt, Direct direct) {
		int offset = direct.getValue();
		if (offset >= 0) {
			return new int[] { pt * dofPerNode + offset };
		} else {
			int[] index = new int[dofPerNode];
			Arrays.setAll(index, i -> pt * dofPerNode + i);
			return index;
		}
	}

	/**
	 * 得到某个单元的自由度
	 * 
	 * @param element 网格单元结点编号
	 * @param direct  方向
	 * @return 单元结点在指定direct方向下的自由度数组
	 */
	public int[] getDof(int[] element, Direct direct) {
		int n = element.length;
		int offset = direct.getValue();
		if (offset >= 0) {
			int[] index = new int[n];
			Arrays.setAll(index, i -> element[i] * dofPerNode + offset);
			return index;
		} else {
			int[] index = new int[n * dofPerNode];
			for (int i = 0; i < n; i++) {
				int start = element[i] * dofPerNode;
				for (int j = 0; j < dofPerNode; j++) {
					index[i * dofPerNode + j] = start + j;
				}
			}
			return index;
		}
	}

	/**
	 * 得到结点集合的自由度
	 * 
	 * @param set    网格结点集合
	 * @param direct 方向
	 * @return 结点集合在指定方向上的自由度
	 */
	int[] getDof(HashSet<Integer> set, Direct direct) {
		return getDof(set, direct, dofPerNode);
	}

	/**
	 * 得到结点集合的自由度
	 * 
	 * @param set    网格结点集合
	 * @param direct 方向
	 * @param dofPerNode 每个结点自由度个数
	 * @return 结点集合在指定方向上的自由度
	 */
	int[] getDof(HashSet<Integer> set, Direct direct, int dofPerNode) {
		int n = set.size();
		int[] index;
		int offset = direct.getValue();
		if (offset >= 0) {
			index = new int[n];
			int i = 0;
			for (int ele : set) {
				index[i++] = ele * dofPerNode + offset;
			}
		} else {
			index = new int[n * dofPerNode];
			int i = 0;
			for (int ele : set) {
				int start = ele * dofPerNode;
				for (int j = 0; j < dofPerNode; j++) {
					index[i * dofPerNode + j] = start + j;
				}
				i++;
			}
		}
		return index;
	}

	/**
	 * 得到单元自由度向量
	 * 
	 * @param ele 单元编号数组
	 * @return 单元自由度向量
	 */
	public int[] getDof(int[] ele) {
		return getDof(ele, dofPerNode);
	}

	/**
	 * 得到单元自由度向量
	 * 
	 * @param ele        单元编号数组
	 * @param dofPerNode 每个结点自由度个数
	 * @return 单元自由度向量
	 */
	public int[] getDof(int[] ele, int dofPerNode) {
		int n = ele.length;
		int[] dof = new int[n * dofPerNode];
		int j, start;
		for (int i = 0; i < n; i++) {
			start = ele[i] * dofPerNode;
			for (j = 0; j < dofPerNode; j++) {
				dof[j + i * dofPerNode] = start + j;
			}
		}
		return dof;
	}

	/**
	 * 由结点集合得到结点对应的整体自由度向量
	 * 
	 * @param ele        结点集合
	 * @param dofPerNode 每个结点自由度个数
	 * @return 结点自由度向量
	 */
	public int[] getDof(HashSet<Integer> ele, int dofPerNode) {
		int n = ele.size();
		ele.toArray();
		int[] dof = new int[n * dofPerNode];
		int j, start;
		int idx = 0;
		for (int myDof : ele) {
			start = myDof * dofPerNode;
			for (j = 0; j < dofPerNode; j++) {
				dof[idx++] = start + j;
			}
		}
		return dof;
	}

	/**
	 * 计算函数在网格结点上的值
	 * 
	 * @param func  标量函数
	 * @param param 标量函数参数
	 * @return 向量, 长度为网格结点个数, 第i个分量表示第i个结点上对应的坐标x下标量函数func的值
	 */
	public double[] valueFromFunc(ScalarFunc func, double[] param) {
		double[] u = new double[mesh.getNv()];
		for (int i = 0, n = u.length; i < n; i++) {
			mesh.getCoordPerNode(i, coordNode);
			u[i] = func.action(coordNode, 0, param);
			// System.out.println(u[i]);
		}
		return u;
	}

	/**
	 * 计算函数在部分网格结点上的值ֵ
	 * 
	 * @param func  标量函数ֵ
	 * @param param 标量函数参数
	 * @param myset 整数结点集合，存储网格结点
	 * @return 向量, 长度为网格结点集合中元素个数, 表示结点集合对应坐标下标量函数func的值ֵ
	 */
	public double[] valueFromSet(ScalarFunc func, double[] param, HashSet<Integer> myset) {
		double[] u = new double[myset.size()];
		int i = 0;
		for (int dof : myset) {
			mesh.getCoordPerNode(dof, coordNode);
			u[i++] = func.action(coordNode, 0, param);
		}
		return u;
	}

	/**
	 * 计算函数在部分网格结点上的值ֵ
	 * 
	 * @param func  标量函数ֵ
	 * @param param 标量函数参数
	 * @param label 网格边界编号
	 * @return 向量, 长度从网格边界编号中提取的网格结点元素个数, 表示这些结点集合对应坐标下标量函数func的值ֵ
	 */
	public double[] valueFromSet(ScalarFunc func, double[] param, int... label) {
		HashSet<Integer> myset = mesh.extractBoundaryNodes(label);
		return valueFromSet(func, param, myset);
	}


/////*********************组装整体矩阵方法(begin)****************************************
	/**
	 * 组装整体矩阵，其中系数以常数的形式给出
	 *
	 * @param constCoef    常数系数数组
	 * @param computeEleMa 单元矩阵函数接口
	 * @param tp           问题类型ֵ
	 * @param A            传入初始矩阵，并返回组装后的矩阵
	 */
	public void assembleGlobalMatrix(double[] constCoef, EleMatFunc computeEleMa, BVPType tp, SRMatrix A, int... domainLabel) {
		double[][] coef = new double[constCoef.length][nPerEle];
		for (int i = 0, n = constCoef.length; i < n; i++) {
			MVO.setValue(coef[i], constCoef[i]);
		}
		if(domainLabel.length == 0) {
			for (int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
				A.assemble(getDofsByElement(i), eleMatrixRegion);
				// System.out.println(Arrays.toString(ele.idx));
			}
			return;
		}
		for(int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			int eleLabel = mesh.getElementLabel(i);
			for(int k = 0, n = domainLabel.length; k < n; k++) {
				if(eleLabel == domainLabel[k]) {
					mesh.getCoordPerElement(ele, coordRegion);
					computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
					A.assemble(getDofsByElement(i), eleMatrixRegion);
					break;
				}
			}
		}
	}

	/**
	 * 组装整体矩阵，其中系数以向量的形式给出
	 *
	 * @param vecCoef      二维数组，行数表示系数个数，列数表示网格结点个数, 在计算单元矩阵时，循环至某个单元(ele_1, ...,
	 *                     ele_N)时， N = nPerElement,
	 *                     从该二维数组的所有行中提取出相应列指标实数vecCoef[i][ele_1], ...,
	 *                     vecCoef[i][ele_N], 将其传入computeEleMa中作为系数计算单元矩阵
	 * @param computeEleMa 单元矩阵函数接口
	 * @param tp           问题类型
	 * @param A            传入初始矩阵，并返回组装后的矩阵
	 */
	public void assembleGlobalMatrix(double[][] vecCoef, EleMatFunc computeEleMa, BVPType tp, SRMatrix A, int... domainLabel) {
		int n = vecCoef.length;
		double[][] coef = new double[n][nPerEle];
		if (domainLabel.length != 0) {
			int lenLabel = domainLabel.length;
			for (int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				int label = mesh.getElementLabel(i);
				for(int j = 0; j < lenLabel; j++) {
					if (label == domainLabel[j]) {
						mesh.getCoordPerElement(ele, coordRegion);
						for (int k = 0; k < n; k++) {
							MVO.getSlice(vecCoef[k], ele, coef[k]);
						}
						computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
						A.assemble(getDofsByElement(i), eleMatrixRegion);
						break;
					}
				}
			}
		} else {
			for (int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				for (int j = 0; j < n; j++) {
					MVO.getSlice(vecCoef[j], ele, coef[j]);
				}
				computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
				A.assemble(getDofsByElement(i), eleMatrixRegion);
			}
		}

	}

	/**
	 * 组装整体矩阵，其中系数以常矩阵形式给出
	 *
	 * @param constMatrixCoef 三维数组, 只用constMatrixCoef[0]中元素作为二维矩阵系数传入单元矩阵函数接口
	 * @param computeEleMa    单元矩阵函数接口
	 * @param tp              问题类型,
	 * @param A               传入初始矩阵，并返回组装后的矩阵
	 */
	public void assembleGlobalMatrix(double[][][] constMatrixCoef, EleMatFunc computeEleMa, BVPType tp, SRMatrix A, int... domainLabel) {
		if (domainLabel.length != 0) {
			int lenLabel = domainLabel.length;
			for (int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				int label = mesh.getElementLabel(i);
				for(int j = 0; j < lenLabel; j++) {
					if (label == domainLabel[j]) {
						mesh.getCoordPerElement(ele, coordRegion);
						computeEleMa.action(ele, coordRegion, constMatrixCoef[0], tp, eleMatrixRegion);
						A.assemble(getDofsByElement(i), eleMatrixRegion);
						break;
					}
				}
			}
		} else {
			for (int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				computeEleMa.action(ele, coordRegion, constMatrixCoef[0], tp, eleMatrixRegion);
				A.assemble(getDofsByElement(i), eleMatrixRegion);
			}
		}
	}

	/**
	 * 组装整体矩阵, 其中系数以标量函数的形式给出
	 * @param coeffFunc
	 * @param param
	 * @param computeEleMa
	 * @param tp
	 * @param A
	 * @param domainLabel
	 */
	public void assembleGlobalMatrix(ScalarFunc[] coeffFunc, double[] param, EleMatFunc computeEleMa, BVPType tp, SRMatrix A, int... domainLabel) {
		double[][] coef = new double[coeffFunc.length][nPerEle];
		if (domainLabel.length != 0) {
			int lenLabel = domainLabel.length;
			for(int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				int eleLabel = mesh.getElementLabel(i);
				for(int j = 0; j < lenLabel; j++) {
					if (eleLabel == domainLabel[j]) {
						mesh.getCoordPerElement(ele, coordRegion);
						getCoefFromFunc(coeffFunc, coordRegion, mesh.getElementLabel(i), param, coef);
						if(MVO.isZero(coef)) {
							continue;
						}
						computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
						A.assemble(getDofsByElement(i), eleMatrixRegion);
						break;
					}
				}
			}
		} else {
			for(int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				getCoefFromFunc(coeffFunc, coordRegion, mesh.getElementLabel(i), param, coef);
				if(MVO.isZero(coef)) {
					continue;
				}
				computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
				A.assemble(getDofsByElement(i), eleMatrixRegion);
			}
		}
	}

	/**
	 * 组装整体矩阵, 其中条数以向量函数的形式给出
	 * @param coeffFunc
	 * @param param
	 * @param computeEleMa
	 * @param tp
	 * @param A
	 * @param domainLabel
	 */
	public void assembleGlobalMatrix(VectorFunc[] coeffFunc, double[][] param, EleMatFunc computeEleMa, BVPType tp, SRMatrix A, int... domainLabel) { }

	/**
	 * 组装整体矩阵, 其中系数以矩阵函数的形式给出
	 * @param matrixCoefFunc
	 * @param param
	 * @param computeEleMa
	 * @param tp
	 * @param A
	 * @param domainLabel
	 */
	public void assembleGlobalMatrix(MatrixFunc matrixCoefFunc, double[][][] param, EleMatFunc computeEleMa, BVPType tp, SRMatrix A, int... domainLabel) {
		if (domainLabel.length == 0) {
			for(int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				double[] ip = mesh.getinteriorPointInElement(coordRegion);
				double[][] coef = matrixCoefFunc.action(ip, mesh.getElementLabel(i), param);
				if(MVO.isZero(coef)) {
					continue;
				}
				computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
				A.assemble(getDofsByElement(i), eleMatrixRegion);
			}
		} else {
			for(int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				int eleLabel = mesh.getElementLabel(i);
				for(int k = 0, n = domainLabel.length; k < n; k++) {
					if(eleLabel == domainLabel[k]) {
						mesh.getCoordPerElement(ele, coordRegion);
						double[] ip = mesh.getinteriorPointInElement(coordRegion);
						double[][] coef = matrixCoefFunc.action(ip, mesh.getElementLabel(i), param);
						if(MVO.isZero(coef)) {
							continue;
						}
						computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
						A.assemble(getDofsByElement(i), eleMatrixRegion);
						break;
					}
				}
			}
		}
	}

	/////////////////////////////////////////////////
	public void assembleStiff(double[] constCoef, BVPType tp, SRMatrix A, int... domainLabel) {
		assembleGlobalMatrix(constCoef, eleStiff, tp, A, domainLabel);
	}
	public void assembleStiff(double[][] vecCoef, BVPType tp, SRMatrix A, int... domainLabel) {
		assembleGlobalMatrix(vecCoef, eleStiff, tp, A, domainLabel);
	}
	public void assembleStiff(double[][][] matCoef, BVPType tp, SRMatrix A, int... domainLabel) {
		assembleGlobalMatrix(matCoef, eleStiff, tp, A, domainLabel);
	}
	public void assembleStiff(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A, int... domainLabel) {
		assembleGlobalMatrix(coefFunc, param, eleStiff, tp, A, domainLabel);
	}
	public void assembleStiff(VectorFunc[] vecFunc, double[][] param, BVPType tp, SRMatrix A, int... domainLabel) {
		assembleGlobalMatrix(vecFunc, param, eleStiff, tp, A, domainLabel);
	}
	public void assembleStiff(MatrixFunc coefFunc, double[][][] param, BVPType tp, SRMatrix A, int... label) {
		assembleGlobalMatrix(coefFunc, param, eleStiff, tp, A, label);
	}
	public void assembleMass(double[] constCoef, BVPType tp, SRMatrix A, int... domainLabel) {
		assembleGlobalMatrix(constCoef, eleMass, tp, A, domainLabel);
	}
	public void assembleMass(double[][] vecCoef, BVPType tp, SRMatrix A, int... domainLabel) {
		assembleGlobalMatrix(vecCoef, eleMass, tp, A, domainLabel);
	}
	public void assembleMass(double[][][] matCoef, BVPType tp, SRMatrix A, int... domainLabel) {
		assembleGlobalMatrix(matCoef, eleMass, tp, A, domainLabel);
	}
	public void assembleMass(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A, int... domainLabel) {
		assembleGlobalMatrix(coefFunc, param, eleMass, tp, A, domainLabel);
	}
	public void assembleMass(VectorFunc[] coefFunc, double[][] param, BVPType tp, SRMatrix A, int... domainLabel) {
		assembleGlobalMatrix(coefFunc, param, eleMass, tp, A, domainLabel);
	}
	public void assembleMass(MatrixFunc coefFunc, double[][][] param, BVPType tp, SRMatrix A, int... domainLabel) {
		assembleGlobalMatrix(coefFunc, param, eleMass, tp, A, domainLabel);
	}
	///////////////////////////////////////////////////
////*********************组装整体矩阵方法(end)****************************************

	/**
	 * 组装整体矩阵，其中系数以标量复值函数的形式给出
	 *
	 * @param coefCFunc    标量复值函数数组，传入相关系数
	 * @param param        标量复值函数参数
	 * @param computeEleMa 计算单元矩阵接口
	 * @param tp           问题类型
	 * @param A            传入初始矩阵，并返回组装后的矩阵
	 */
	public void assembleGlobalMatrix(ScalarCFunc[] coefCFunc, Complex[] param, EleMatCFunc computeEleMa, BVPType tp, SCMatrix A) {
		Complex[][] coef = Complex.zeros(coefCFunc.length, nPerEle);
		Complex[][] eleMatrixRegion = Complex.zeros(dofPerElement, dofPerElement);
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			getCoefFromFunc(coefCFunc, coordRegion, mesh.getElementLabel(i), param, coef);
			computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
			A.assemble(getDof(ele), eleMatrixRegion);
		}
	}

	/**
	 *
	 * @param coefFunc
	 * @param param
	 * @param eleMat
	 * @param tp
	 * @param A
	 * @param domainLabel
	 */
	public void assembleGlobalMatrix(ScalarFunc[] coefFunc, double[] param, EleMMatFunc eleMat, BVPType[] tp, SRMatrix[] A, int... domainLabel) {
		double[][] coef = new double[coefFunc.length][nPerEle];
		double[][][] eleMatrix = new double[A.length][dofPerElement][dofPerElement];
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			getCoefFromFunc(coefFunc, coordRegion, mesh.getElementLabel(i), param, coef);
			eleMat.action(ele, coordRegion, coef, tp, eleMatrix);
			int[] dof = getDofsByElement(i);
			for (int j = 0, n = A.length; j < n; j++) {
				A[j].assemble(dof, eleMatrix[j]);
			}
		}
	}

	public void assembleGlobalMatrix(MatrixCFunc matrixCoefFunc, double[][][] param, EleMatCFunc computeEleMa, BVPType tp,
									 SCMatrix A, int... domainLabel) {
		//double[][] coef = new double[coeffFunc.length][nPerEle];
		Complex[][] eleMatrixRegion = Complex.zeros(dofPerElement, dofPerElement);
		for(int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			int eleLabel = mesh.getElementLabel(i);
			for(int k = 0, n = domainLabel.length; k < n; k++) {
				if(eleLabel == domainLabel[k]) {
					mesh.getCoordPerElement(ele, coordRegion);
					double[] ip = mesh.getinteriorPointInElement(coordRegion);
					Complex[][] coef = matrixCoefFunc.action(ip, mesh.getElementLabel(i), param);
					computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
					A.assemble(getDofFromDofs(ele), eleMatrixRegion);
				}
			}
		}
	}

	/**
	 * 组装整体矩阵, 参数以数组形式给出, 表示给出的
	 * @param constCoef
	 * @param computeEleMa
	 * @param tp
	 * @param A
	 * @param domainLabel
	 */
	public void assembleGlobalMatrix(double[] constCoef, EleMatFunc computeEleMa, BVPType tp, SCMatrix A, int... domainLabel) {
		double[][] coef = new double[constCoef.length][nPerEle];
		for (int i = 0, n = constCoef.length; i < n; i++) {
			MVO.setValue(coef[i], constCoef[i]);
		}
		if(domainLabel == null) {
			for (int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
				A.assemble(getDofsByElement(i), eleMatrixRegion);
			}
		} else {
			for(int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				int eleLabel = mesh.getElementLabel(i);
				for(int k = 0, n = domainLabel.length; k < n; k++) {
					if(eleLabel == domainLabel[k]) {
						mesh.getCoordPerElement(ele, coordRegion);
						computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
						A.assemble(getDofsByElement(i), eleMatrixRegion);
						break;
					}
				}
			}
		}
	}
	
	public void assembleGlobalMatrix(Complex[] constCoef, EleMatCFunc computeEleMa, BVPType tp, SCMatrix A) {
		Complex[][] coef = Complex.zeros(constCoef.length, nPerEle);
		Complex[][] eleMatrixRegion = Complex.zeros(dofPerElement, dofPerElement);
		for (int i = 0, n = constCoef.length; i < n; i++) {
			for (int j = 0; j < nPerEle; j++) {
				coef[i][j].real = constCoef[i].real;
				coef[i][j].image = constCoef[i].image;
			}
		}
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
			A.assemble(getDofFromDofs(ele), eleMatrixRegion);
		}
	}

	/**
	 * 组装整体矩阵，其中系数以向量的形式给出
	 * 
	 * @param vecCoef      二维数组，行数表示系数个数，列数表示网格结点个数, 在计算单元矩阵时，循环至某个单元(ele_1, ..., ele_N)时， N = nPerElement,
	 *                                                       从该二维数组的所有行中提取出相应列指标实数vecCoef[i][ele_1], ...,
	 *                     vecCoef[i][ele_N], 将其传入computeEleMa中作为系数计算单元矩阵
	 * @param computeEleMa 单元矩阵函数接口
	 * @param tp           问题类型
	 * @param A            传入初始矩阵，并返回组装后的矩阵
	 */
	public void assembleGlobalMatrix(Complex[][] vecCoef, EleMatCFunc computeEleMa, BVPType tp, SCMatrix A) {
		int n = vecCoef.length;
		Complex[][] coef = Complex.zeros(n, nPerEle);
		Complex[][] eleMatrixRegion = Complex.zeros(dofPerElement, dofPerElement);
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			for (int j = 0; j < n; j++) {
				Complex.getSlice(vecCoef[j], ele, coef[j]);
			}
			computeEleMa.action(ele, coordRegion, coef, tp, eleMatrixRegion);
			A.assemble(getDof(ele), eleMatrixRegion);
		}
	}


/////********************组装整体向量(begin)********************
	/**
	 * 组装整体向量, 参数以标量的形式给出
	 *
	 * @param constCoef    标量参数数组
	 * @param computeEleVe 单元向量
	 * @param tp           问题类型
	 * @param RHS          组装的右端项
	 */
	public void assembleGlobalVector(double[] constCoef, EleVecFunc computeEleVe, BVPType tp, double[] RHS, int... domainLabel) {
		double[][] coef = new double[constCoef.length][nPerEle];
		for (int i = 0, n = constCoef.length; i < n; i++) {
			MVO.setValue(coef[i], constCoef[i]);
		}
		int lenLabel = domainLabel.length;
		if (lenLabel != 0) {
			for (int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				int label = mesh.getElementLabel(i);
				for(int j = 0; j < lenLabel; j++) {
					if (label == domainLabel[j]) {
						mesh.getCoordPerElement(ele, coordRegion);
						computeEleVe.action(ele, coordRegion, coef, tp, eleVecRegion);
						assembleVector(eleVecRegion, getDof(ele), RHS);
						break;
					}
				}
			}
		} else {
			for (int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				computeEleVe.action(ele, coordRegion, coef, tp, eleVecRegion);
				assembleVector(eleVecRegion, getDof(ele), RHS);
			}
		}
	}

	public void assembleGlobalVector(double[][] vecCoef, EleVecFunc computeEleVe, BVPType tp, double[] RHS, int... domainLabel) {
		double[][] coef = new double[vecCoef.length][nPerEle];
		int lenLabel = domainLabel.length;
		if (lenLabel != 0) {
			for (int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				int label = mesh.getElementLabel(i);
				for(int j = 0; j < lenLabel; j++) {
					if (label == domainLabel[j]) {
						mesh.getCoordPerElement(ele, coordRegion);
						for (int jj = 0, n1 = vecCoef.length; jj < n1; jj++) {
							for (int kk = 0; kk < nPerEle; kk++) {
								coef[jj][kk] = vecCoef[jj][ele[kk]];
							}
						}
						if(MVO.isZero(coef)) {
							continue;
						}
						computeEleVe.action(ele, coordRegion, coef, tp, eleVecRegion);
						assembleVector(eleVecRegion, getDof(ele), RHS);
						break;
					}
				}

			}
		} else {
			for (int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				for (int j = 0, n1 = vecCoef.length; j < n1; j++) {
					for (int k = 0; k < nPerEle; k++) {
						coef[j][k] = vecCoef[j][ele[k]];
					}
				}
				computeEleVe.action(ele, coordRegion, coef, tp, eleVecRegion);
				assembleVector(eleVecRegion, getDof(ele), RHS);
			}
		}
	}

	/**
	 * 组装整体向量，其中系数以标量函数的形式给出
	 *
	 * @param coeffFunc    标量函数数组，传入相关系数
	 * @param param        标量函数参数
	 * @param computeEleVe 单元向量函数接口
	 * @param tp           问题类型
	 * @param RHS          传入初始向量，并返回组装后的向量
	 */
	public void assembleGlobalVector(ScalarFunc[] coeffFunc, double[] param, EleVecFunc computeEleVe, BVPType tp, double[] RHS, int... domainLabel) {
		double[][] coef = new double[coeffFunc.length][nPerEle];
		int lenLabel = domainLabel.length;
		if(lenLabel != 0) {
			for (int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				int label = mesh.getElementLabel(i);
				for(int j = 0; j < lenLabel; j++) {
					if(label == domainLabel[j]) {
						mesh.getCoordPerElement(ele, coordRegion);
						getCoefFromFunc(coeffFunc, coordRegion, label, param, coef);
						if (MVO.isZero(coef)) {
							continue;
						}
						computeEleVe.action(ele, coordRegion, coef, tp, eleVecRegion);
						assembleVector(eleVecRegion, getDof(ele), RHS);
						break;
					}
				}
			}
		} else {
			for (int i = 0; i < nE; i++) {
				// elementInfo.get(i, ele, label, coord);
				// mesh.getInfoPerElement(i, eleRegion, label, coordRegion);
				int[] ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				getCoefFromFunc(coeffFunc, coordRegion, mesh.getElementLabel(i), param, coef);
				computeEleVe.action(ele, coordRegion, coef, tp, eleVecRegion);
				assembleVector(eleVecRegion, getDof(ele), RHS);
			}
		}
	}

	/**
	 * 组装整体向量, 其中系数以标量形式给出, 并传入已知向量参数, 可用于组装热载荷向量， 其中系数为热模量,
	 * 以常数的形式给出，已知向量参数为已经计算得到的温度向量
	 *
	 * @param constCoef 标量函数数组，传入相关系数
	 * @param coef1    向量参数, 以矩阵形式给出, 行数为向量个数, 列数为网格结点个数
	 * @param eleFunc  计算单元向量接口数组
	 * @param tp       问题类型
	 * @param RHS      传入初始向量，并返回组装后的向量
	 */
	public void assembleGlobalVector(double[] constCoef, double[][] coef1, EleVecFunc eleFunc, BVPType tp, double[] RHS, int... domainLabel) {
		int len1 = constCoef.length;
		int len2 = coef1.length;
		double[][] coef = new double[len1 + len2][nPerEle];
		for (int i = 0, n = constCoef.length; i < n; i++) {
			MVO.setValue(coef[i], constCoef[i]);
		}
		int lenLabel = domainLabel.length;
		if (lenLabel == 0) {
			for (int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				for (int j = 0; j < len2; j++) {
					MVO.getSlice(coef1[j], ele, coef[len1 + j]);
				}
				eleFunc.action(ele, coordRegion, coef, tp, eleVecRegion);
				assembleVector(eleVecRegion, getDofsByElement(i), RHS);
			}
		} else {
			for (int i = 0; i < nE; i++) {
				int[] ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				int label = mesh.getElementLabel(i);
				for(int j = 0; j < lenLabel; j++) {
					if (label == domainLabel[j]) {
						for (int k = 0; j < len2; j++) {
							MVO.getSlice(coef1[k], ele, coef[len1 + j]);
						}
						eleFunc.action(ele, coordRegion, coef, tp, eleVecRegion);
						assembleVector(eleVecRegion, getDofsByElement(i), RHS);
						break;
					}
				}
			}
		}
	}

	/**
	 * 组装整体向量，其中系数以标量函数的形式给出
	 * 
	 * @param coeffFunc    标量函数数组，传入相关系数
	 * @param param        标量函数参数
	 * @param computeEleVe 单元向量函数接口
	 * @param tp           问题类型
	 * @param RHS          传入初始向量，并返回组装后的向量
	 */
	public void assembleGlobalVector(ScalarCFunc[] coeffFunc, Complex[] param, EleVecCFunc computeEleVe, BVPType tp, Complex[] RHS) {
		Complex[][] coef = Complex.zeros(coeffFunc.length, nPerEle);
		Complex[] eleVecRegion = Complex.zeros(dofPerElement);
		for (int i = 0; i < nE; i++) {
			// elementInfo.get(i, ele, label, coord);
			// mesh.getInfoPerElement(i, eleRegion, label, coordRegion);
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			getCoefFromFunc(coeffFunc, coordRegion, mesh.getElementLabel(i), param, coef);
			computeEleVe.action(ele, coordRegion, coef, tp, eleVecRegion);
			assembleCVector(eleVecRegion, getDof(ele), RHS);
		}
	}

	/**
	 * 组装整体向量，其中系数以标量函数的形式给出
	 * 
	 * @param coeffFunc     标量函数数组，传入相关系数
	 * @param param        标量函数参数
	 * @param computeEleVe 计算单元向量接口数组
	 * @param tp           问题类型数组, 与eleMat数组长度相同，标识每一个单元向量函数接口的问题类型
	 * @param RHS            初始二维数组, 行数为eleMat数组长度, 分别将各单元向量加入到相应的二维数组各行中
	 */
	public void assembleGlobalVector(ScalarFunc[] coeffFunc, double[] param, EleMVecFunc computeEleVe, BVPType tp,
			double[][] RHS) {
		double[][] coef = new double[coeffFunc.length][nPerEle];
		int len = RHS.length;
		eleMultiVecRegion = new double[len][dofPerElement];
		for (int i = 0; i < nE; i++) {
			// elementInfo.get(i, ele, label, coord);
			// mesh.getInfoPerElement(i, eleRegion, label, coordRegion);
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			getCoefFromFunc(coeffFunc, coordRegion, mesh.getElementLabel(i), param, coef);
			computeEleVe.action(ele, coordRegion, coef, tp, eleMultiVecRegion);
			int[] dof = getDof(ele);
			for (int j = 0; j < len; j++) {
				assembleVector(eleMultiVecRegion[j], dof, RHS[j]);
			}
		}
	}

	/**
	 * 组装整体向量, 其中系数以标量形式给出, 并传入已知向量参数, 可用于组装热载荷向量， 其中系数为热模量,
	 * 以函数的形式给出，已知向量参数为已经计算得到的温度向量
	 * 
	 * @param coefFunc 标量函数数组，传入相关系数
	 * @param param    标量函数参数
	 * @param coef1    向量参数, 以矩阵形式给出, 行数为向量个数, 列数为网格结点个数
	 * @param eleFunc  计算单元向量接口数组
	 * @param tp       问题类型
	 * @param RHS      传入初始向量，并返回组装后的向量
	 */
	public void assembleGlobalVector(ScalarFunc[] coefFunc, double[] param, double[][] coef1, EleVecFunc eleFunc,
			BVPType tp, double[] RHS) {
		int len1 = coefFunc.length;
		int len2 = coef1.length;
		double[][] coef = new double[len1 + len2][nPerEle];
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			getCoefFromFunc(coefFunc, coordRegion, mesh.getElementLabel(i), param, coef);
			for (int j = 0; j < len2; j++) {
				MVO.getSlice(coef1[j], ele, coef[len1 + j]);
			}
			eleFunc.action(ele, coordRegion, coef, tp, eleVecRegion);
			assembleVector(eleVecRegion, getDof(ele), RHS);
		}
	}

	/**
	 * 组装整体向量, 其中系数以标量形式给出, 并传入两个已知向量参数,
	 * 
	 * @param coeffFunc
	 * @param param
	 * @param coef1
	 * @param coef2
	 * @param eleFunc
	 * @param tp
	 * @param RHS
	 */
	public void assembleGlobalVector(ScalarFunc[] coeffFunc, double[] param, double[][] coef1, double[][] coef2,
			EleMMVecFunc eleFunc, BVPType tp, double[][] RHS) {
		double[][] coef = new double[coeffFunc.length][nPerEle];
		int len = RHS.length;
		eleMultiVecRegion = new double[len][dofPerElement];
		int n = coef1.length;
		double[][] ccoef1 = new double[n][dofPerElement];
		for (int i = 0; i < nE; i++) {
			// elementInfo.get(i, ele, label, coord);
			// mesh.getInfoPerElement(i, eleRegion, label, coordRegion);
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			getCoefFromFunc(coeffFunc, coordRegion, mesh.getElementLabel(i), param, coef);
			int[] dof = getDof(ele);
			for (int j = 0; j < n; j++) {
				MVO.getSlice(coef1[j], dof, ccoef1[j]);
			}
			eleFunc.action(ele, coordRegion, coef, ccoef1, coef2, tp, eleMultiVecRegion);
			for (int j = 0; j < len; j++) {
				assembleVector(eleMultiVecRegion[j], dof, RHS[j]);
			}
		}
	}

	public void assembleGlobalVector(MatrixFunc coeffFunc, double[][][] param, double[][] coef1, double[][] coef2,
			EleMMVecFunc eleFunc, BVPType tp, double[][] RHS) {
		// double[][] coef = new double[coeffFunc.length][nPerEle];
		int len = RHS.length;
		eleMultiVecRegion = new double[len][dofPerElement];
		int n = coef1.length;
		double[][] ccoef1 = new double[n][dofPerElement];
		for (int i = 0; i < nE; i++) {
			// elementInfo.get(i, ele, label, coord);
			// mesh.getInfoPerElement(i, eleRegion, label, coordRegion);
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			// getCoefFromFunc(coeffFunc, coordRegion, ele.label, param, coef);
			double[] ip = mesh.getinteriorPointInElement(coordRegion);
			double[][] coef = coeffFunc.action(ip, mesh.getElementLabel(i), param);
			int[] dof = getDof(ele);
			for (int j = 0; j < n; j++) {
				MVO.getSlice(coef1[j], dof, ccoef1[j]);
			}
			eleFunc.action(ele, coordRegion, coef, ccoef1, coef2, tp, eleMultiVecRegion);
			// int[] dof = getDof(ele.idx);
			for (int j = 0; j < len; j++) {
				assembleVector(eleMultiVecRegion[j], dof, RHS[j]);
			}
		}
	}

	public void assembleGlobalVector(MatrixFunc coeffFunc, ScalarFunc[] sfunc, double[][][] param, double[][] coef1,
			double[][] coef2, EleMMMVecFunc eleFunc, BVPType tp, double[][] RHS) {
		// double[][] coef = new double[coeffFunc.length][nPerEle];
		int len = RHS.length;
		eleMultiVecRegion = new double[len][dofPerElement];
		int sLen = sfunc.length;
		double[][] sCoef = new double[sLen][dofPerElement];
		int n = coef1.length;
		double[][] ccoef1 = new double[n][dofPerElement];
		double[] param2 = new double[1];
		for (int i = 0; i < nE; i++) {
			// elementInfo.get(i, ele, label, coord);
			// mesh.getInfoPerElement(i, eleRegion, label, coordRegion);
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			int label = mesh.getElementLabel(i);
			double[] ip = mesh.getinteriorPointInElement(coordRegion);
			double[][] coef = coeffFunc.action(ip, label, param);
			getCoefFromFunc(sfunc, coordRegion, label, param2, sCoef);
			int[] dof = getDof(ele);
			for (int j = 0; j < n; j++) {
				MVO.getSlice(coef1[j], dof, ccoef1[j]);
			}
			eleFunc.action(ele, coordRegion, coef, sCoef, ccoef1, coef2, tp, eleMultiVecRegion);
			// int[] dof = getDof(ele.idx);
			for (int j = 0; j < len; j++) {
				assembleVector(eleMultiVecRegion[j], dof, RHS[j]);
			}
		}
	}

	public void assembleGlobalVector(Complex[] constCoef, EleVecCFunc computeEleVe, BVPType tp, Complex[] RHS) {
		Complex[][] coef = Complex.zeros(constCoef.length, nPerEle);
		for (int i = 0, n = constCoef.length; i < n; i++) {
			for (int j = 0; j < nPerEle; j++) {
				coef[i][j].real = constCoef[i].real;
				coef[i][j].image = constCoef[i].image;
			}
		}
		Complex[] eleVecRegion = new Complex[dofPerElement];
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			// mesh.getInfoPerElement(i, eleRegion, label, coordRegion);
			computeEleVe.action(ele, coordRegion, coef, tp, eleVecRegion);
			assembleCVector(eleVecRegion, getDof(ele), RHS);
		}
	}
	

	public void assembleGlobalVector(MatrixFunc coeffFunc, double[][][] param, EleVecFunc computeEleVe, BVPType tp, double[] RHS) {
		// double[][] coef = new double[coeffFunc.length][nPerEle];
		for (int i = 0; i < nE; i++) {
			// elementInfo.get(i, ele, label, coord);
			// mesh.getInfoPerElement(i, eleRegion, label, coordRegion);
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			// getCoefFromFunc(coeffFunc, coordRegion, ele.label, param, coef);
			double[] ip = mesh.getinteriorPointInElement(coordRegion);
			double[][] coef = coeffFunc.action(ip, mesh.getElementLabel(i), param);
			computeEleVe.action(ele, coordRegion, coef, tp, eleVecRegion);
			assembleVector(eleVecRegion, getDof(ele), RHS);
		}
	}

	public void assembleGlobalVector(MatrixFunc coeffFunc, double[][][] param, EleMVecFunc computeEleVe, BVPType tp, double[][] RHS) {
		int len = RHS.length;
		eleMultiVecRegion = new double[len][dofPerElement];
		for (int i = 0; i < nE; i++) {
			// elementInfo.get(i, ele, label, coord);
			// mesh.getInfoPerElement(i, eleRegion, label, coordRegion);
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			double[] ip = mesh.getinteriorPointInElement(coordRegion);
			double[][] coef = coeffFunc.action(ip, mesh.getElementLabel(i), param);
			computeEleVe.action(ele, coordRegion, coef, tp, eleMultiVecRegion);
			int[] dof = getDof(ele);
			for (int j = 0; j < len; j++) {
				assembleVector(eleMultiVecRegion[j], dof, RHS[j]);
			}
		}
	}
////********************组装整体向量(end)********************

	public void assembleBoundaryMatrix(ScalarFunc[] coefFunc, double[] param, EleMatFunc computeBdMatrix, BVPType tp,
			SRMatrix A, int... bdLabel) {
		double[][] coef = new double[coefFunc.length][nPerBoundary];
		for (int i = 0; i < nB; i++) {
			int[] ele = mesh.getBoundary(i);
			mesh.getCoordPerBoundary(ele, coordBoundary);
			int label = mesh.getBoundaryLabel(i);
			// mesh.getInfoPerBoundary(i, eleBoundary, label, coordBoundary);
			for (int j = 0, nLabel = bdLabel.length; j < nLabel; j++) {
				if (label == bdLabel[j]) {
					getCoefFromFunc(coefFunc, coordBoundary, label, param, coef);
					computeBdMatrix.action(ele, coordBoundary, coef, tp, eleMatrixBoundary);
					A.assemble(getDof(ele), eleMatrixBoundary);
					break;
				}
			}
		}
	}

	public void assembleBoundaryMatrix(ScalarCFunc[] coefFunc, Complex[] param, EleMatCFunc computeBdMatrix, BVPType tp,
			SCMatrix A, int... bdLabel) {
		Complex[][] coef = Complex.zeros(coefFunc.length, nPerBoundary);
		Complex[][] eleMatrixBoundary = Complex.zeros(nPerBoundary, dofPerBoundary);
		for (int i = 0; i < nB; i++) {
			int[] ele = mesh.getBoundary(i);
			mesh.getCoordPerBoundary(ele, coordBoundary);
			int label = mesh.getBoundaryLabel(i);
			// mesh.getInfoPerBoundary(i, eleBoundary, label, coordBoundary);
			for (int j = 0, nLabel = bdLabel.length; j < nLabel; j++) {
				if (label == bdLabel[j]) {
					getCoefFromFunc(coefFunc, coordBoundary, label, param, coef);
					computeBdMatrix.action(ele, coordBoundary, coef, tp, eleMatrixBoundary);
					A.assemble(getDof(ele), eleMatrixBoundary);
					break;
				}
			}
		}
	}

	public void assembleBoundaryMatrix(double[] constCoef, EleMatFunc computeBdMatrix, BVPType tp, SRMatrix A,
			int... bdLabel) {
		int n = constCoef.length;
		double[][] coef = new double[n][nPerBoundary];
		for (int i = 0; i < n; i++) {
			MVO.setValue(coef[i], constCoef[i]);
		}
		for (int i = 0; i < nB; i++) {
			int[] ele = mesh.getBoundary(i);
			mesh.getCoordPerBoundary(ele, coordBoundary);
			int label = mesh.getBoundaryLabel(i);
			for (int j = 0, nLabel = bdLabel.length; j < nLabel; j++) {
				if (label == bdLabel[j]) {
					computeBdMatrix.action(ele, coordBoundary, coef, tp, eleMatrixBoundary);
					A.assemble(getDof(ele), eleMatrixBoundary);
					break;
				}
			}
		}
	}
	
	public void assembleBoundaryMatrix(Complex[] constCoef, EleMatCFunc computeBdMatrix, BVPType tp, SCMatrix A,
			int... bdLabel) {
		int n = constCoef.length;
		Complex[][] coef = Complex.zeros(n, nPerBoundary);
		Complex[][] eleMatrixBoundary = Complex.zeros(dofPerBoundary, dofPerBoundary);
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < nPerEle; j++) {
				coef[i][j].real = constCoef[i].real;
				coef[i][j].image = constCoef[i].image;
			}
		}
		for (int i = 0; i < nB; i++) {
			int[] ele = mesh.getBoundary(i);
			mesh.getCoordPerBoundary(ele, coordBoundary);
			int label = mesh.getBoundaryLabel(i);
			for (int j = 0, nLabel = bdLabel.length; j < nLabel; j++) {
				if (label == bdLabel[j]) {
					computeBdMatrix.action(ele, coordBoundary, coef, tp, eleMatrixBoundary);
					A.assemble(getDof(ele), eleMatrixBoundary);
					break;
				}
			}
		}
	}

	public void assembleBoundaryMatrix(double[][] vecCoef, EleMatFunc computeBdMatrix, BVPType tp, SRMatrix A, int... bdLabel) {
		int n = vecCoef.length;
		double[][] coef = new double[n][nPerBoundary];
		for (int i = 0; i < nB; i++) {
			int[] ele = mesh.getBoundary(i);
			mesh.getCoordPerBoundary(ele, coordBoundary);
			for (int j = 0; j < n; j++) {
				MVO.getSlice(vecCoef[j], ele, coef[j]);
			}
			int label = mesh.getBoundaryLabel(i);
			for (int j = 0, nLabel = bdLabel.length; j < nLabel; j++) {
				if (label == bdLabel[j]) {
					computeBdMatrix.action(ele, coordBoundary, coef, tp, eleMatrixBoundary);
					A.assemble(getDof(ele), eleMatrixBoundary);
					break;
				}
			}
		}
	}

	/**
	 * ͨ���ṩϵ��������װ�߽�����
	 * 
	 * @param coeffFunc    ϵ����������
	 * @param param        ϵ����������
	 * @param computeEleVe ��Ԫ��������
	 * @param tp           ��Ԫ��������
	 * @param RHS          �����Ҷ�����
	 * @param bdLabel      �߽�������
	 */
	public void assembleBoundaryVector(ScalarFunc[] coeffFunc, double[] param, EleVecFunc computeEleVe, BVPType tp,
			double[] RHS, int... bdLabel) {
		double[][] coef = new double[coeffFunc.length][nPerBoundary];
		for (int i = 0; i < nB; i++) {
			int[] ele = mesh.getBoundary(i);
			mesh.getCoordPerBoundary(ele, coordBoundary);
			// mesh.getInfoPerBoundary(i, eleBoundary, label, coordBoundary);
			int label = mesh.getBoundaryLabel(i);
			for (int j = 0, nLabel = bdLabel.length; j < nLabel; j++) {
				if (label == bdLabel[j]) {
					getCoefFromFunc(coeffFunc, coordBoundary, label, param, coef);
					computeEleVe.action(ele, coordBoundary, coef, tp, eleVecBoundary);
					assembleVector(eleVecBoundary, getDof(ele), RHS);
					break;
				}
			}
		}
	}
	
	public void assembleBoundaryVector(ScalarCFunc[] coeffFunc, Complex[] param, EleVecCFunc computeEleVe, BVPType tp,
			Complex[] RHS, int... bdLabel) {
		Complex[][] coef = Complex.zeros(coeffFunc.length, nPerBoundary);
		Complex[] eleVecBoundary = Complex.zeros(dofPerBoundary);
		for (int i = 0; i < nB; i++) {
			int[] ele = mesh.getBoundary(i);
			mesh.getCoordPerBoundary(ele, coordBoundary);
			// mesh.getInfoPerBoundary(i, eleBoundary, label, coordBoundary);
			int label = mesh.getBoundaryLabel(i);
			for (int j = 0, nLabel = bdLabel.length; j < nLabel; j++) {
				if (label == bdLabel[j]) {
					getCoefFromFunc(coeffFunc, coordBoundary, label, param, coef);
					computeEleVe.action(ele, coordBoundary, coef, tp, eleVecBoundary);
					assembleCVector(eleVecBoundary, getDof(ele), RHS);
					break;
				}
			}
		}
	}
	

	public void assembleBoundaryVector(double[] constCoef, EleVecFunc computeEleVe, BVPType tp, double[] RHS,
			int... bdLabel) {
		double[][] coef = new double[constCoef.length][nPerBoundary];
		for (int i = 0, n = constCoef.length; i < n; i++) {
			MVO.setValue(coef[i], constCoef[i]);
		}
		for (int i = 0; i < nB; i++) {
			int[] ele = mesh.getBoundary(i);
			mesh.getCoordPerBoundary(ele, coordBoundary);
			int label = mesh.getBoundaryLabel(i);
			// mesh.getInfoPerBoundary(i, eleBoundary, label, coordBoundary);
			for (int j = 0, nLabel = bdLabel.length; j < nLabel; j++) {
				if (label == bdLabel[j]) {
					// getCoefFromFunc(coeffFunc, coordBoundary, ele.label, param, coef);
					computeEleVe.action(ele, coordBoundary, coef, tp, eleVecBoundary);
					// System.out.println("OK");
					assembleVector(eleVecBoundary, getDof(ele), RHS);
					break;
				}
			}
		}
	}
	
	public void assembleBoundaryVector(Complex[] constCoef, EleVecCFunc computeEleVe, BVPType tp, Complex[] RHS,
			int... bdLabel) {
		Complex[][] coef = Complex.zeros(constCoef.length, nPerBoundary);
		for (int i = 0, n = constCoef.length; i < n; i++) {
			for (int j = 0; j < nPerEle; j++) {
				coef[i][j].real = constCoef[i].real;
				coef[i][j].image = constCoef[i].image;
			}
		}
		Complex[] eleVecBoundary = Complex.zeros((nPerBoundary));
		for (int i = 0; i < nB; i++) {
			int[] ele = mesh.getBoundary(i);
			mesh.getCoordPerBoundary(ele, coordBoundary);
			int label = mesh.getBoundaryLabel(i);
			// mesh.getInfoPerBoundary(i, eleBoundary, label, coordBoundary);
			for (int j = 0, nLabel = bdLabel.length; j < nLabel; j++) {
				if (label == bdLabel[j]) {
					// getCoefFromFunc(coeffFunc, coordBoundary, ele.label, param, coef);
					computeEleVe.action(ele, coordBoundary, coef, tp, eleVecBoundary);
					// System.out.println("OK");
					assembleCVector(eleVecBoundary, getDof(ele), RHS);
					break;
				}
			}
		}
	}
	

	public void assembleBoundaryVector(double[][] vecCoef, EleVecFunc computeEleVe, BVPType tp, double[] RHS,
			int... bdLabel) {
		double[][] coef = new double[vecCoef.length][nPerBoundary];
		for (int i = 0; i < nB; i++) {
			int[] ele = mesh.getBoundary(i);
			mesh.getCoordPerBoundary(ele, coordBoundary);
			int label = mesh.getBoundaryLabel(i);
			// mesh.getInfoPerBoundary(i, eleBoundary, label, coordBoundary);
			for (int j = 0, nLabel = bdLabel.length; j < nLabel; j++) {
				if (label == bdLabel[j]) {
					// getCoefFromFunc(coeffFunc, coordBoundary, ele.label, param, coef);
					for (int ki = 0, n = vecCoef.length; ki < n; ki++) {
						for (int kj = 0; kj < nPerBoundary; kj++) {
							coef[ki][kj] = vecCoef[ki][ele[kj]];
						}
					}
					computeEleVe.action(ele, coordBoundary, coef, tp, eleVecBoundary);
					assembleVector(eleVecBoundary, getDof(ele), RHS);
					break;
				}
			}
		}
	}

	public void map(double[] source, FES fesDestination, double[] destination) {

		ParamCheck.checkEqual(fesDestination.GetNdof(), destination.length);
	}

	public void assembleHeatThetaScheme0(ScalarFunc[] coefFunc, double[] param, BVPType[] tp, EleMatFunc elementStiff,
			EleMatFunc elementMass, SRMatrix LMatrix, SRMatrix RMatrix) {
		double[][] kxx = new double[1][nPerEle];
		double[][] rho = new double[1][nPerEle];
		double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
		double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
		double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
		double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
		double dt = param[0];
		double theta = param[1];
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			int label = mesh.getElementLabel(i);
			for (int j = 0, nj = nPerEle; j < nj; j++) {
				kxx[0][j] = coefFunc[0].action(coordRegion[j], label, param);
				rho[0][j] = coefFunc[1].action(coordRegion[j], label, param);
			}
			elementStiff.action(ele, coordRegion, kxx, tp[0], eleMatrixStiff);
			elementMass.action(ele, coordRegion, rho, tp[1], eleMatrixMass);
			for (int ii = 0; ii < dofPerElement; ii++) {
				for (int jj = 0; jj < dofPerElement; jj++) {
					eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
					eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
				}
			}
			LMatrix.assemble(ele, eleLMatrix);
			RMatrix.assemble(ele, eleRMatrix);
		}
	}

	public void assembleHeatThetaScheme0(double[] coefFunc, BVPType[] tp, EleMatFunc elementStiff,
			EleMatFunc elementMass, SRMatrix LMatrix, SRMatrix RMatrix) {
		double[][] kxx = new double[1][nPerEle];
		double[][] rho = new double[1][nPerEle];
		double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
		double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
		double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
		double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
		Arrays.fill(kxx[0], coefFunc[0]);
		Arrays.fill(rho[0], coefFunc[1]);
		double dt = coefFunc[2];
		double theta = coefFunc[3];
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			elementStiff.action(ele, coordRegion, kxx, tp[0], eleMatrixStiff);
			elementMass.action(ele, coordRegion, rho, tp[1], eleMatrixMass);
//			System.out.println(MVO.toString(eleMatrixMass));
//			System.out.println(MVO.toString(eleMatrixStiff));
//			System.exit(0);
			for (int ii = 0; ii < dofPerElement; ii++) {
				for (int jj = 0; jj < dofPerElement; jj++) {
					eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
					eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
				}
			}
			LMatrix.assemble(ele, eleLMatrix);
			RMatrix.assemble(ele, eleRMatrix);
		}
	}


	public double[] computeReaction(ScalarFunc[] coeffFunc, double[] param, BVPType tp, EleMatFunc eleFunc, double[] solution) {
		double[][] coef = new double[coeffFunc.length][nPerEle];
		double[] solDof = new double[dofPerElement];
		double[] result = new double[solution.length];
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			int[] dof = getDof(ele);
			mesh.getCoordPerElement(ele, coordRegion);
			getCoefFromFunc(coeffFunc, coordRegion, mesh.getElementLabel(i), param, coef);
			eleFunc.action(ele, coordRegion, coef, tp, eleMatrixRegion);
			MVO.getSlice(solution, dof, solDof);
			eleFunc.action(ele, coordRegion, coef, tp, eleMatrixRegion);
			double[] reaction = MVO.matmul(eleMatrixRegion, solDof);
			assembleVector(reaction, dof, result);
			// System.out.println(Arrays.toString(ele.idx));
		}
		return result;
	}

	public double[] computeReaction(double[] constCoef, BVPType tp, EleMatFunc eleFunc, double[] solution) {
		double[][] coef = new double[constCoef.length][nPerEle];
		double[] solDof = new double[dofPerElement];
		double[] result = new double[solution.length];
		for (int i = 0, n = constCoef.length; i < n; i++) {
			MVO.setValue(coef[i], constCoef[i]);
		}
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			int[] dof = getDof(ele);
			mesh.getCoordPerElement(ele, coordRegion);
			eleFunc.action(ele, coordRegion, coef, tp, eleMatrixRegion);
			MVO.getSlice(solution, dof, solDof);
			eleFunc.action(ele, coordRegion, coef, tp, eleMatrixRegion);
			double[] reaction = MVO.matmul(eleMatrixRegion, solDof);
			assembleVector(reaction, dof, result);
			// System.out.println(Arrays.toString(ele.idx));
		}
		return result;
	}

	public double[] computeReaction(double[][] vecCoef, BVPType tp, EleMatFunc eleFunc, double[] solution) {
		int n = vecCoef.length;
		double[][] coef = new double[n][nPerEle];
		double[] solDof = new double[dofPerElement];
		double[] result = new double[solution.length];
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			int[] dof = getDof(ele);
			mesh.getCoordPerElement(ele, coordRegion);
			for (int j = 0; j < n; j++) {
				MVO.getSlice(vecCoef[j], ele, coef[j]);
			}
			MVO.getSlice(solution, dof, solDof);
			eleFunc.action(ele, coordRegion, coef, tp, eleMatrixRegion);
			double[] reaction = MVO.matmul(eleMatrixRegion, solDof);
			assembleVector(reaction, dof, result);
		}
		return result;
	}


	public void assembleSource(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS) {
		assembleGlobalVector(coefFunc, param, eleSource, tp, RHS);
	}

	public void assembleSource(double[] constCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(constCoef, eleSource, tp, RHS);
	}

	public void assembleSource(double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(vecCoef, eleSource, tp, RHS);
	}

	public void assembleSource(double[][] vecCoef, BVPType tp, double[] RHS, int... label) {
		assembleGlobalVector(vecCoef, eleSource, tp, RHS, label);
	}

	public void assembleFlux(ScalarFunc[] coeffFunc, double[] param, BVPType tp, double[] RHS, int... bdLabel) {
		assembleBoundaryVector(coeffFunc, param, eleFlux, tp, RHS, bdLabel);
	}

	public void assembleFlux(double[] constCoef, BVPType tp, double[] RHS, int... bdLabel) {
		assembleBoundaryVector(constCoef, eleFlux, tp, RHS, bdLabel);
	}

	public void assembleFlux(double[][] vecCoef, BVPType tp, double[] RHS, int... bdLabel) {
		assembleBoundaryVector(vecCoef, eleFlux, tp, RHS, bdLabel);
	}

	/**
	 * 有限元空间关联的网格
	 */
	public Mesh mesh;

	/**
	 * 自由度个数
	 */
	protected int ndof;

	/**
	 * 问题维数，等于网格维数
	 */
	protected int nDim;

	/**
	 * 每个单元结点数
	 */
	protected int nPerEle;

	/**
	 * 每个边界结点数
	 */
	protected int nPerBoundary;

	/**
	 * 每个结点自由度个数
	 */
	protected int dofPerNode;

	/**
	 * 单元个数
	 */
	protected int nE;

	/**
	 * 边界个数
	 */
	protected int nB;

	/**
	 * 每个单元自由度数
	 */
	protected int dofPerElement;

	/**
	 * 每个边界自由度数
	 */
	protected int dofPerBoundary;

	/**
	 * 形成的稀疏矩阵半带宽
	 */
	protected int nBand;

	/**
	 * 结点坐标数组
	 */
	double[] coordNode;

	/**
	 * 单元坐标数组
	 */
	double[][] coordRegion;

	/**
	 * 边界单元数组
	 */
	double[][] coordBoundary;

	/**
	 * 单元矩阵数组
	 */
	double[][] eleMatrixRegion;

	/**
	 * 边界单元矩阵数组
	 */
	double[][] eleMatrixBoundary;

	/**
	 * 单元向量数组
	 */
	double[] eleVecRegion;

	double[][] eleMultiVecRegion;

	/**
	 * 边界向量数组
	 */
	double[] eleVecBoundary;

	/**
	 * 结点自由度数组
	 */
	int[][] dofs;

	EleMatFunc eleStiff;
	EleMatFunc eleMass;
	EleVecFunc eleSource;
	EleVecFunc eleFlux;
	final int[] DIRECT_X = {0};
	final int[] DIRECT_Y = {1};
	final int[] DIRECT_Z = {2};
	final int[] DIRECT_XY = {0, 1};
	final int[] DIRECT_YZ = {1, 2};
	final int[] DIRECT_XZ = {0, 2};
	final int[] DIRECT_XYZ = {0, 1, 2};
}
